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•  A  hybrid  strategy  is  proposed  for  mechanistic  modeling  of  fuel  cell  stacks. 

•  The  strategy  is  verified  with  the  full  set  of  equations  for  a  10-cell  stack. 

•  The  strategy  preserves  geometrical  resolution  and  captures  essential  physics. 

•  In  particular,  captures  redistribution  of  current  and  heat  flux  across  the  cells. 

•  80%  reductions  in  computational  cost;  allows  for  simulation  of  large  stacks. 
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In  general,  detailed  mechanistic  models  for  fuel  cell  stacks  that  seek  to  capture  the  local  transport 
phenomena  are  computationally  expensive.  In  this  context,  we  propose  a  hybrid  modeling  strategy,  in 
which  the  steady-state  conservation  equations  are  solved  iteratively  in  two  separate  groups:  The  first 
comprises  asymptotically-reduced  governing  equations  for  momentum,  mass,  and  species,  which  are 
solved  as  a  transient-like  propagation  problem;  the  second  comprises  the  full  set  of  equations  for  energy 
and  charge,  which  are  solved  as  an  elliptic  stationary  problem.  Physically,  the  segregation  is  justified  by 
the  nature  of  the  dependent  variables;  in  essence,  the  first  group  covers  local  variables  on  the  cell  level 
and  the  second  involves  global  variables  on  the  stack  level.  We  demonstrate  the  methodology  for  a 
steady-state  detailed  mechanistic  model  of  a  proton  exchange  membrane  fuel  cell  stack  comprising  2  to 
350  cells  subjected  to  non-uniform  operating  conditions  across  the  cells;  e.g.,  a  stack  comprising  350 
cells  takes  less  than  an  hour  to  solve.  The  proposed  methodology  is  generic  and  can  also  be  employed  for 
other  systems  where  transport  phenomena  occur  on  different  length  scales  and  involve  some  form  of 
slenderness. 

©  2013  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Mathematical  modeling  and  simulation  have  found  widespread 
use  in  the  research  and  development  of  fuel  cells;  however,  the 
majority  of  detailed  mechanistic  models  that  provide  geometrical 
resolution  and  resolve  the  essential  physics  have  focused  on  the  cell 
level  so  far  [1-8],  At  the  stack  level,  only  a  few  detailed  mechanistic 
models  can  be  found  [9-18],  which  are  limited  to  small  stacks  of  up 
to  around  5  to  10  cells.  Modeling  of  larger  stacks  generally  involves 
simplifications  [19—37]  with  a  loss  in  the  level  of  detail  and  reso¬ 
lution  of  the  salient  features  of  the  electrochemical  and  transport 
phenomena. 
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In  essence,  detailed  mechanistic  models  typically  consider 
conservation  of  mass,  momentum,  species,  energy  and  charge  in 
the  form  of  a  system  of  coupled  nonlinear  elliptic  partial  differential 
equations  (PDEs)  with  the  relevant  boundary  conditions  and 
constitutive  relations  to  ensure  a  well-posed  problem.  While  these 
models  can  be  solved  reasonably  fast  and  efficiently  for  a  single  cell, 
the  associated  computational  cost  for  stacks  quickly  becomes 
prohibitive  because  of  the  large  number  of  functional  domains- 
separator  plates  (sp),  gas  flow  fields  (ff),  coolant  flow  fields  (cff), 
porous  backings  (pb),  catalyst  layer  (cl),  electrolyte  (el)  —  in  each 
and  every  cell  in  a  stack  (see  Fig.  1 ).  One  strategy  to  deal  with  such 
large  simulations  is  to  deploy  clusters  of  computers  or  other  types 
of  parallel  computing  environments  to  speed  up  calculations,  yet 
these  simulations  still  take  on  the  order  of  hours  or  days  for  fuel  cell 
stacks  [1,7,12,25,28,38—40],  Another  possible  strategy  to  reduce  the 
overall  complexity  and  associated  computational  cost  of  a  stack 
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Fig.  1.  Schematic  for  a  PEMFC  stack  comprising  n  cells,  denoted  by  j  (a);  mathematical  nature  of  the  governing  equations  for  the  full  stack  model  (b)  and  reduced  stack  model  (c): 
elliptic  PDEs  (H),  parabolic  PDEs  (  — )  and  ODEs  (|). 


model  is  to  address  the  underlying  mathematical  nature  of  the 
system  of  equations.  Care  has  to  be  taken  though  so  as  not  to  reduce 
the  fidelity  of  model  predictions  especially  when  reductions  sac¬ 
rifice  leading-order  phenomena  or  dimensionality.  This  strategy 
was  demonstrated  in  Refs.  [41],  where  a  leading-order  reduced 
mathematical  model  for  a  proton  exchange  membrane  fuel  cell 
(PEMFC)  stack  was  derived  and  shown  to  cut  down  the  computa¬ 
tional  cost  by  2-3  orders  of  magnitude  (depending  on  stack  size). 
The  leading-order  solution  was,  however,  found  to  be  unable  to 
capture  the  redistribution  of  physical  quantities  for  stacks  sub¬ 
jected  to  significant  perturbations  between  cells— a  key  phenome¬ 
non  in  fuel  cell  stacks  [42-44], 

In  light  of  the  promise  the  reduced  model  [41  ]  holds  in  terms  of 
computational  efficiency  even  for  large  stacks  of  hundreds  of  cells 
or  more,  we  shall  here  seek  to  overcome  its  main  limitation  by 
exploring  a  hybrid  strategy  that  exploits  the  inherent  nature  of  the 
dependent  variables.  In  short,  we  shall  argue  that  the  dependent 
variables  exist  either  locally  in  each  cell  or  globally  throughout  the 
stack  whence  their  PDEs  can  be  split  up  into  two  separate  groups. 
The  local  variables  can  now  be  solved  at  leading  order  with  a  space¬ 
marching  algorithm  whilst  the  global  are  allowed  to  retain  their 
elliptic  nature;  iteration  between  the  two  groups  then  links  them 
together  into  a  complete  stack  model.  We  demonstrate  the  pro¬ 
posed  strategy  for  the  steady-state  operation  of  a  10-cell  PEMFC 
stack  subject  to  non-uniform  operating  conditions  for  the  indi¬ 
vidual  cells. 

The  layout  of  the  paper  is  as  follows.  First,  we  introduce  the 
mathematical  formulation  for  the  full  and  proposed  hybrid  model 
for  a  PEMFC  stack,  which  is  then  followed  by  a  summary  of  the 
hybrid  coupling  methodology.  After  a  description  of  numerics 
necessary  to  solve  the  models,  the  results  of  the  hybrid  model  for 
both  the  global  and  local  behavior  are  verified  against  those  of  the 
full  model.  The  savings  in  computational  cost  in  terms  of  memory 
usage  and  computing  time  are  then  discussed.  Finally,  we  draw 
conclusions  and  highlight  how  the  present  approach  can  be 
adapted  to  other  systems. 


2.  Mathematical  formulation 

Let  us  start  by  turning  our  attention  towards  a  typical  fuel  cell 
stack  as  illustrated  schematically  in  Fig.  1,  which  comprises  several 
single  cells  connected  in  series;  each  cell  further  contains  several 
functional  layers.  A  detailed  mechanistic  model  needs  to  consider 
coupled  transport  phenomena  —  mass,  momentum,  species,  en¬ 
ergy,  and  charge  -  in  all  of  the  relevant  length  scales  that  can  be 
found  in  the  stack.  Since  the  cells  in  a  stack  are  connected  via 
impermeable  graphite  or  metallic  plates,  transport  of  mass,  mo¬ 
mentum,  and  species  is  confined  to  individual  cells;  in  contrast,  the 
transport  of  energy  and  charge  occurs  across  the  cells  and 
throughout  the  entire  stack.  Thus,  the  transport  phenomena  in  the 
fuel  cell  stack,  from  the  viewpoint  of  their  scales  of  occurrence,  can 
be  classified  into  two  categories: 

•  Local  to  a  cell:  conservation  of  mass,  momentum  and  species.  At 
this  scale,  the  inherent  slenderness  in  the  functional  layers  of 
the  cell  and  the  relatively  impermeable  nature  of  all  layers 
except  for  the  flow  field  can  be  exploited  to  reduce  the  elliptic 
PDEs  to  parabolic  counterparts  and  even  ordinary  differential 
equations  (ODEs)  with  significant  savings  in  computational 
cost  [41,45-48], 

•  Global  to  the  stack:  conservation  of  energy  and  charge.  At  this 
scale,  we  retain  the  elliptic  nature  of  the  PDEs  since  there  is  no 
slenderness  that  can  be  exploited  at  leading-order.  (N.B.:  It 
should  be  possible  to  still  reduce  the  equations,  but  we  shall 
not  pursue  that  here;  see  the  Conclusions  section  for  a  short 
discussion.) 

This  in  turn  suggests  that  a  hybrid  modeling  strategy  based  on 
these  two  groups  should  be  feasible  and  allow  for  significant 
reduction  in  computational  cost  without  sacrificing  the  fidelity  of 
model  predictions  for  the  stack.  We  shall  justify  the  strategy  a 
posteriori  by  comparing  the  model  predictions  from  the  hybrid 
strategy  with  solutions  from  the  full  non-reduced  set  of  equations. 
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To  demonstrate  the  proposed  hybrid  methodology,  we  consider 
a  two-dimensional  (2D)  PEMFC  stack  system  comprising  n  number 
of  single  cells  and  n  -  1  number  of  liquid  coolant  plates  (Fig.  1 );  and 
the  single-phase  steady-state  model  [41,45]  that  accounts  for 
conservation  of  mass,  momentum,  species  (O2,  H2O  and  N2  in  the 
cathode;  H2,  H2O  and  N2  in  the  anode),  energy,  and  charge  (see 
Appendix  A).  The  cells  are  equipped  with  porous-type  flow  fields 
for  reactant  gases  as  well  as  for  coolant  plates  and  are  operating  in 
coflow  mode.  We  use  porous  flow  fields  in  order  to  exploit  the  fact 
that  they  can  conveniently  be  reduced  from  three  dimensions  (3D) 
to  two  dimensions  (2D)  and  note  that  the  other  types  of  common 
flow  fields  -  parallel  and  serpentine  flow  fields  -  can  also  be 
represented  by  porous  counterparts  through  spatial  smoothing 
[25,48,49],  As  discussed  earlier,  the  governing  equations  for  con¬ 
servation  of  mass,  momentum,  and  species  can  be  reduced  to  a  set 
of  parabolic  PDEs  in  the  flow  fields,  and  ODEs  in  the  other  func¬ 
tional  layers  (see  Appendix  B).  We  refer  to  [45]  for  detailed  dis¬ 
cussion  on  the  reduction  methodology  and  validation  with 
experiments;  and  [41,45]  for  other  model  details  -  boundary 
conditions,  constitutive  equations  and  base  case  parameters.  (N.B.: 
The  strategy  is  generic  and  can  be  extended  to  account  for  multi¬ 
phase  flow  or  3D  geometries  that  cannot  be  reduced  to  2D.) 

3.  Hybrid  coupling  methodology 

To  couple  the  solutions  for  the  dependent  variables  at  the  local 
and  global  level,  we  iterate  between  the  two  as  illustrated  in  the 
flowchart  in  Fig.  2.  In  essence,  the  following  steps  are  carried  out: 

1  The  local  dependent  variables  in  the  form  of  reduced  equations 
are  solved;  here,  we  solve  reduced  equations  for  the  global 
variables  as  well  to  find  good  initial  guesses  for  the  elliptic 
solver  (not  necessary,  but  does  speed  up  overall  calculations). 


Fig.  2.  Flowchart  for  hybrid  simulation  strategy  for  fuel  cell  stacks. 


2  The  distributions  of  the  local  dependent  variables  are  mapped 
to  a  geometry  and  computational  mesh  for  the  solver  of  the 
global  dependent  variables:  energy  and  charge  transport. 

3  The  global  dependent  variables  are  solved. 

4  The  global  dependent  variables  are  mapped  to  the  geometry 
and  computational  mesh  for  the  local  dependent  variables. 

5  Step  1  is  repeated  without  solving  for  temperature  and  charge, 
which  are  evaluated  from  the  mapped  solution  obtained  by 
step  4. 

6  Step  2-5  are  repeated  until  the  relative  error  between  suc¬ 
cessive  iterations  for  each  variable  is  less  than  the  specified 
tolerance. 

7  Finally,  the  solutions  for  all  dependent  variables  are  mapped  to 
geometry  for  postprocessing. 


4.  Numerics 

The  commercial  finite-element  solver,  Comsol  3.5a  [50],  which 
allows  for  the  solution  of  generic  differential  equations,  was 
employed  to  solve  the  2D  full  and  hybrid  models.  The  elliptic 
equations  in  the  full  model  as  well  as  in  the  hybrid  model  were 
solved  with  an  elliptic  solver  that  resolves  both  the  space  di¬ 
mensions  altogether,  whereas  the  reduced  set  of  equations  (para¬ 
bolic  PDEs  and  ODEs)  was  solved  by  resolving  they-direction  as  the 
computational  domain  and  then  marching  along  the  x-direction 
with  a  space  marching  algorithm  using  a  transient  solver.  In  other 
words,  the  local  dependent  variables  were  solved  on  a  ID  transient 
domain  and  the  global  variables  on  a  2D  stationary  domain. 

A  total  of  62  mesh  elements  were  used  in  the  y-direction  for 
each  cell  in  the  stack  (both  in  full  and  reduced  model)  and  60  el¬ 
ements  in  the  x-direction  for  the  full  set,  whereas  the  reduced 
model  is  simulated  for  120  time  steps  in  the  x-direction  treating  x 
as  a  timelike  variable.  The  direct  solver,  Pardiso  was  chosen  as  the 
linear  solver  for  the  full  non-reduced  equations  in  both  the  full  and 
hybrid  model.  For  the  reduced  set  of  equations,  we  chose  Umfpack  as 
the  linear  solver  and  backward  differentiation  formula  with  adap¬ 
tive  time  stepping  as  the  space  marching  scheme.  Quadratic 
Lagrange  elements  were  employed  with  a  relative  convergence 
tolerance  of  10  4  for  all  variables.  Further,  solutions  from  all  the 
models  were  tested  for  mesh  independence. 

We  employed  automated  model  generation  [41]  to  build  the 
numerical  stack  models  both  for  full  and  reduced  set  of  equations 
by  exploiting  the  bidirectional  interface  between  Comsol  and  Matlab 
[51  ].  In  brief,  the  automated  procedure  in  the  form  of  a  Matlab  script 
manipulates  CoMSOL-associated  structures,  called  the  FEM  struc¬ 
tures  that  contain  the  numerical  formulation  of  the  entire  model, 
and  automates  the  steps  that  are  necessary  to  obtain  numerical 
solutions  of  the  mathematical  model:  viz.,  drawing  the  geometry, 
meshing,  implementing  the  mathematical  equations,  solving,  and 
post-processing.  This  approach  thus  allows  for  a  large  reduction  in 
time,  since  the  time  spent  on  manually  setting  up  and  solving  large 
stacks  is  shortened  considerably. 

The  hybrid  stack  model  couples  the  two  FEM  structures  ac¬ 
cording  to  the  methodology  described  in  the  last  section.  For 
bijective  mapping  of  the  dependent  variables  from  one  domain  to 
another,  the  same  number  of  elements  (mesh  points)  are  assigned 
in  the  y-direction,  and  the  output  times  for  the  transient  solver  are 
chosen  accordingly.  Further,  the  solutions  on  the  Lagrange  nodes 
are  stored  in  different  order  for  the  two  domains:  as  a  vector  for 
2D  stationary  and  as  a  matrix  for  ID  transient.  Hence,  we  permute 
the  solution  from  the  source  domain  to  the  order  of  the  destina¬ 
tion  domain  while  mapping  the  solutions.  The  convergence  cri¬ 
terion  for  the  hybrid  strategy  was  computed  as  the  weighted 
Euclidean  norm  [50], 
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~£(N/Wf)2, 


(i) 


where  N  is  the  number  of  degrees  of  freedom,  Ej  is  the  estimated 
absolute  error  between  iterations,  Wj  is  defined  as  max((Jj,Sj), 
where  U-,  is  the  current  solution,  and  Si  was  set  to  0.1  of  the  mean  of 
the  solution. 

All  computations  were  carried  out  on  a  workstation  with  a 
hexacore  processor  of  12  threads  (3.2  GHz)  and  a  total  of  64  GB 
random  access  memory  (RAM).  The  peak  memory  usage  was  esti¬ 
mated  from  the  task  manager  taking  the  difference  of  idle  com¬ 
puter  memory  and  the  peak  memory  during  the  simulation  runs. 
The  real  execution  times  (wall-clock  time)  were  determined  from 
Matlab’s  tic-toc  commands  with  all  unnecessary  processes  stopped 
to  ensure  reasonably  accurate  estimates. 

Finally,  we  note  that  the  hybrid  strategy  does  not  rely  on  the 
finite-element  as  such;  the  strategy  could  be  implemented  with 
other  numerical  techniques  as  well,  such  as  the  finite-difference 
and  finite-volume  method.  Furthermore,  one  could  implement 
the  reduced  and  full  set  of  equations  in  a  high-level  programming 
language  such  as  Fortran  or  C;  and  implement  numerical  algo¬ 
rithms  based  on,  e.g.,  a  Keller  box  or  a  modified  box  method. 

5.  Verification 

So  far,  we  have  employed  an  understanding  of  transport  phe¬ 
nomena  scales  —  local  to  a  cell  or  global  to  the  stack  —  to  propose  a 
hybrid-modeling  and  simulation  strategy  for  fuel  cell  stacks.  Now, 
in  order  to  verify  that  the  hybrid  model  does  indeed  capture  the 
behavior  of  the  full  model  with  redistribution  of  current  and  heat 
flux  for  the  stack,  we  proceed  by  considering  a  10-cell  PEMFC  stack 
operating  under  non-identical  operating  conditions.  In  this  regard, 
we  introduce  a  leading-order  perturbation  in  the  inlet  velocities  for 
the  cathodes,  given  by 

U§  =  U*n  +  jAU,  (2) 

where  j  denotes  the  number  of  the  cell  in  the  stack,  U'cnj  is  the  inlet 
velocity  at  the  cathode  of  cell  j,  L/jn  =  1  ms-1,  and  AU  is  the  incre¬ 
ment  in  inlet  velocity  from  cell  to  cell.  For  demonstration  purposes, 


Fig.  4.  Local  current  density  distributions  for  a  10-cell  stack  (£stack  =  3  V)  along  the  x- 
axis  at  the  interface  between  the  cathode  catalyst  layer  and  the  membrane  for  the  full 
set  of  equations  in  cell  ( ▲ )  1,  (•)  5,  ( ■ )  10;  and  corresponding  predictions  of  the  hybrid 


AU  is  chosen  such  that  the  inlet  velocity  increases  by  a  factor  of  2 
between  the  first  and  the  last  cell,  which  is  larger  than  the  typical 
variation  in  inlet  velocities  between  5  and  25%  that  have  been 
studied  for  stacks  comprising  25  to  100  cells  [21,22,52], 

We  start  by  comparing  the  global  behavior  in  the  form  of  po¬ 
larization  curves  predicted  by  the  models  based  on  the  full  and 
hybrid  sets  of  equations,  as  shown  in  Fig.  3.  Overall,  the  agreement 
is  good  with  a  maximum  relative  error  less  than  1%.  Similarly,  and 
most  importantly,  the  redistribution  of  the  current  density  and  the 
temperature  distribution  between  cells  in  the  stack  due  to  the 
perturbed  inlet  velocities  at  the  cathodes  is  captured  with  a 
maximum  relative  error  less  than  1%;  see  Figs.  4  and  5  for  the  first, 
fifth,  and  tenth  cell  in  the  stack  at  an  operating  voltage  of  3  V.  The 


Fig.  5.  Local  temperature  distributions  for  a  10-cell  stack  (£stack  =  3  V)  along  the  x-axis 
at  the  interface  between  the  cathode  catalyst  layer  and  the  membrane  for  the  full  set  of 
equations  in  cell  ( A )  1,  (•}  5,  ( ■ )  10;  and  corresponding  predictions  of  the  hybrid  set  in 
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Fig.  6.  Local  oxygen  concentration  distributions  for  a  10-cell  stack  (Estack  =  3  V)  along 
the  x-axis  at  the  interface  between  the  cathode  catalyst  layer  and  the  membrane  for 
the  full  set  of  equations  in  cell  ( A )  1,  (•)  5,  ( ■ )  10;  and  corresponding  predictions  of 
the  hybrid  set  in  lines. 
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Number  of  cells 

Fig.  7.  Scale-up  tests  for  the  hybrid  model:  (■)  memory  usage,  (•)  time  for  setting  up 
the  numerical  stack  models,  and  ( A )  convergence  time  (at  a  typical  operating  voltage 
of  roughly  0.6  V  for  each  cell). 


stack  voltage  is  chosen  to  be  in  the  concentration  polarization  re¬ 
gion,  which  is  the  most  nonlinear  operating  regime  that  corre¬ 
sponds  to  high  current  densities  with  larger  spatial  variations  in  the 
dependent  variables.  The  hybrid  strategy  has  thus  removed  the 
limitation  of  the  earlier  derived  stack  model  [41],  which  could  not 
capture  the  redistribution  across  the  cells. 

Finally,  let  us  verily  the  local  dependent  variables  that  were 
solved  with  reduced  governing  equations.  We  do  so  for  the  oxygen 
concentration  in  Fig.  6  and  again  find  good  agreement  with  the 
predicted  distribution  from  the  full  set  of  non-reduced  equations; 
the  other  dependent  variables  show  the  same  agreement,  but  are 
not  shown  here  for  the  sake  of  brevity. 

6.  Computational  cost  and  efficiency 

Now  that  the  fidelity  of  the  hybrid  model  has  been  verified,  let 
us  proceed  to  address  the  computational  cost  and  scalability.  We 
expect  that  the  hybrid  strategy  should  lead  to  significant  compu¬ 
tational  savings  as  compared  to  simply  solving  the  full  set  of 
equations  both  in  terms  of  memory  and  convergence  time.  This  is 
indeed  the  case,  as  can  be  seen  in  Table  1,  which  compares  the 
degrees  of  freedom  (DoF),  memory  requirement  and  the  time 
required  to  reach  the  converged  solution  for  the  full  model  and 
hybrid  counterpart  for  three  different  stacks,  comprising  2,  5,  and 
10  cells,  respectively.  All  stacks  are  simulated  for  the  same  amount 


of  perturbation  in  cathode  inlet  velocities  across  the  stack,  i.e.,  the 
inlet  velocity  increases  by  a  factor  of  2  between  the  first  and  the  last 
cell.  The  DoFs  for  the  hybrid  strategy  decrease  as  compared  to  the 
full  set  of  equations  which  leads  to  substantial  savings  in  the 
convergence  time  and  memory  requirements.  Generally  speaking, 
we  find  that  around  80%  less  memory  and  execution  time  is  needed 
for  hybrid  simulations  than  for  full  model  simulations  for  stack 
sizes  up  to  around  10  cells:  e.g.,  a  10-cell  stack  at  6  V  can  be 
simulated  within  1  min  and  2  GB  of  RAM  with  the  hybrid  strategy 
as  compared  to  6  min  and  9  GB  of  RAM  for  the  full  counterpart.  As 
the  operating  potential  decreases,  the  mathematical  equations 
become  more  nonlinear  rendering  the  model  more  difficult  to 
solve,  which  is  evident  from  the  increased  computational  time  for 
both  the  models  at  lower  voltages. 

On  our  workstation,  the  computational  cost  for  the  full  model 
becomes  prohibitive  for  stacks  with  50  or  more  cells,  unless  one 
switches  to  iterative  solvers  that  are  less  robust  and  need  fine- 
tuning.  However,  the  hybrid  strategy  can  simulate  the  behavior  of 
stacks  up  to  around  350  cells  with  a  robust  solver  as  depicted  in 
Fig.  7.  Here,  we  see  that  the  computational  cost  in  terms  of  memory 
usage,  time  required  for  setting  up  the  models  (done  automatically 
with  the  automated  code  generator),  and  convergence  time  in¬ 
crease  roughly  linearly  with  the  number  of  cells.  One  could  pursue 
larger  stacks  by  switching  to  an  iterative  solver  for  the  global  var¬ 
iables  governed  by  elliptic  PDEs  in  order  to  decrease  the  memory 
requirements— doing  so  would  be  at  the  expense  of  robustness, 
however. 


Computational  cost  for  the  full  and  hybrid  sets;  the  numbers  in  the  brackets  indicate 
the  time  required  to  automatically  generate  the  numerical  stack  model  before 
solving  it. 


10-Cell  stack 


DoF 

Time  (s)  for  1.6, 1.2,  0.8  V 
Memory  (GB) 

DoF 

Time  (s)  for  4, 3,  2  V 
Memory  (GB) 

DoF 

Memory  (GB) 


1.5  x  105 
(3)  39,57, 110 
2.2 

3.9  x  105 

(6)  130,150,  340 

5 

7.8  x  105 
(11)210,330,  750 


6.1  x  104 
(3)  9,13,18 
0.67 

(7)  26,27,39 


j.l  A  1U 

(13)  35,53,  82 


7.  Conclusions 

We  have  proposed  and  demonstrated  a  hybrid  strategy  to  tackle 
the  prohibitive  computational  costs  of  solving  detailed  mechanistic 
models  for  fuel  cell  stacks  whilst  capturing  the  essential  physics-in 
particular  the  redistribution  of  current  and  heat  flux  between  cells 
in  a  perturbed  stack. 

In  essence,  the  strategy  relies  on  exploiting  the  scales  and  nature 
of  the  dependent  variables:  chemical  species,  momentum  and  mass 
are  local  to  cells  and  their  governing  equations  can  be  reduced 
mathematically  by  invoking  the  slenderness  and  impermeable 
nature  of  the  functional  layers  found  at  the  cell  level;  whereas 
charge  and  temperature  are  global  to  the  stack,  and  thus  retain 
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their  elliptic  nature.  The  strategy  was  verified  with  good  agreement 
by  comparing  with  solutions  obtained  from  solving  the  non- 
reduced  full  set  of  equations. 

Around  80%  computational  savings  both  in  terms  of  memory 
usage  and  execution  time  were  achieved  with  the  hybrid  strategy 
for  stacks  up  to  around  10  cells.  In  addition,  the  strategy  allows  for 
the  simulation  of  larger  stacks:  e.g.,  it  took  less  than  an  hour  to 
simulate  a  perturbed  350-cell  stack.  The  scalability  and  associated 
low  computational  cost  should  allow  for  various  applied  studies  on 
fuel  cell  stacks  e.g.,  modeling  of  statistical  and/or  random  variations 
and  perturbations  in  the  operating  conditions  and  material  prop¬ 
erties  arising  during  manufacture  and  assembly,  establishing 
quality  control  guidelines  for  rejection  or  acceptance  of  a  cell  dur¬ 
ing  assembly,  and  optimizing  operating  and  design  parameters  for  a 
stack. 

Here,  we  solve  the  full  set  of  equations  for  energy  and  charge 
transport  coupled  with  the  reduced  equations  for  other  transport 
processes.  Another  possible  approach  could  be  to  continue  with  the 
asymptotically  reduced  equations  for  all  the  transport  phenomena 
augmented  with  the  next-order  solutions  of  their  asymptotic  series 
to  resolve  the  redistribution  of  the  energy  and  charge  between  the 
cells.  An  alternate  strategy  could  be  to  secure  analytical  solutions  of 
the  elliptic  governing  equations  for  charge  and  energy  and  couple 
these  with  the  reduced  formulation  for  the  local  dependent 
variables. 

On  a  final  note,  the  classification  of  transport  phenomena  based 
on  their  scales  and  associated  model  reductions  for  the  local  vari¬ 
ables  is  a  fundamental  concept  and  could  thus  be  extended  to  3D 
models  and  other  electrochemical  systems,  such  as  batteries,  flow 
batteries,  and  supercapacitors.  Furthermore,  the  strategy  is  not 
limited  to  electrochemical  stacks,  but  could  also  be  applied  to  other 
systems  that  exhibit  some  degree  of  slenderness  for  local  variables: 
e.g.,  in  catalytic  monolith  reactors  comprising  multiple  slender 
channels  where  the  mass,  momentum,  and  species  transport  are 
limited  to  single  channels  but  the  heat  transport  occurs  throughout 
the  reactor. 
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Appendix  A.  Full  set  of  governing  equations 

We  consider  conservation  of  mass,  momentum,  species,  energy 


and  charge  in  the  PEMFC,  expressed  as 

v-(pv)  =  Smass,  (ff,  pb,  cl,  el)  (A.l) 

Vp  = -^v  +  Smom,(ff,  pb,  cl,  el)  (A.2) 

V-Nj  =  Si,  (ff,  pb,  cl,  el)  (A.3) 

pCpV-VT  =  V'(/ceffVT)  +Sten  ip,  (everywhere)  (A.4) 

V-im  =  Spot, (cl,  el)  (A.5) 

V  is  =  -Spot,(sp,  ff,  pb,  cl)  (A.6) 


where  the  molar  fluxes  of  species,  N,  and  current  densities,  i  are 
defined  as 

No2  =  Co2 v  —  Dq^Vco2  ,  (ff,  pb,  cl  in  the  cathode)  (A.7) 

nh2  =  Ch2v  -  D^Vch2,  (ff,  pb,  cl  in  the  anode)  (A.8) 

^-fcDH2o,mvA  (el)  (Ag) 

Ch2oV  -  £>HfoVcH2o  (elsewhere) 

im  =  (A.10) 

is  =  — ersV0s-  (A.11) 

For  the  coolant  flow  fields,  we  only  need  to  consider  conserva¬ 
tion  of  energy  and  charge,  which  can  be  written  as 

p(l)Cp)Ucoolex- VT  =  V-(keffVT)  +  Stamp,  (cff)  (A.12) 

V-(— <rsV0s)  =  0  (cff)  (A.13) 

In  the  above  equations,  p  is  the  density,  v  is  the  velocity,  p  is  the 
pressure,  p  is  the  dynamic  viscosity,  k  is  the  hydraulic  permeability 
of  the  porous  medium,  cp  is  the  specific  heat  capacity,  T  is  the 
temperature,  kefris  the  effective  thermal  conductivity,  im  is  the  ionic 
current  density,  and  is  is  the  electronic  current  density.  In  the  molar 
flux  of  species  (Eqs.  A.7  and  A.8),  Ci  and  denote  the  molar 
concentration  and  effective  diffusivity  of  species  i.  The  flux  of  water 
in  the  membrane  (Eq.  (A.9))  due  to  electro-osmotic  drag  and 
diffusion  is  expressed  using  a  phenomenological  model  [53]  in 
terms  of  the  membrane  water  content,  X  where  rip  is  electroosmotic 
drag  coefficient,  im  is  the  current  density  carried  by  protons,  F  is 
Faraday’s  constant,  DH20m  is  the  diffusivity  of  water  in  the  mem¬ 
brane,  pm  and  Mm  are  the  density  and  equivalent  weight  of  the  dry 
membrane,  respectively.  In  the  flux  of  current  density  (Eqs.  A.10 
and  A.11),  <pm  represents  the  potential  of  the  ionic  phase  and  rps 
the  solid  phase;  om  and  as  are  the  electrical  conductivities  of  proton 
and  electron  transport,  respectively.  In  Eq.  (A.12),  pm  and  c®  are  the 
density  and  the  specific  heat  capacity  of  the  liquid  coolant  (H2O  in 
this  case),  If 001  is  the  average  velocity  of  the  flow  through  the 
porous  coolant  plates  (N.B.  for  passive  flow  in  a  channel  comprising 
a  porous  medium  with  slip-conditions,  a  constant  velocity  profile  is 
obtained,  which  is  given  by  If001  in  Eq.  (A.12)),  and  ex  is  the 
streamwise  coordinate  vector. 


Appendix  B.  Reduced  equations  for  conservation  of  mass, 
momentum,  and  species 

We  solve  for  the  conservation  of  mass  and  momentum,  given  by 


dp  p  9^ 

_  <W9VA2  (ff) 

(B.1) 

0x  ~  Kp  9y 

<"> 

(B.2) 

=  5™ 

s,  (pb,  cl,  el) 

(B.3) 
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9 P 

ey 


(B.4) 


The  streamfunctions,  \j/,  are  defined  through 


u 


lcty 

p0y‘ 


(B.5) 


such  that  they  satisfy  conservation  of  mass  in  the  flow  fields;  here, 
u  and  v  are  the  velocities  in  the  x-  and  y-direction,  respectively. 
Note  that  we  solve  for  i p  and  p  in  the  flow  fields,  in  lieu  of  u,  v,  and  p, 
whereas  we  solve  for  v  and  p  in  the  gas  diffusion  layer,  catalyst 
layers  and  membranes.  The  second  term  in  right  hand  side  of  Eq.  B.l 
originates  from  the  momentum  source,  Smom;  refer  to  [41]  for 
details. 

The  reduced  equations  for  the  conservation  of  species  can  be 


expressed  as 

|(C,U)  +  A(q,)  -  (ff) 

(B-6) 

<Pb,d) 

(B-7) 

(B.8) 

List  of  symbols 


cF 

form-drag  constant 

Ci 

molar  concentration  of  species  i,  mol  m 

Cp 

specific  heat  capacity,  J  kg-1  K-1 

Di 

diffusivity  of  species  i,  m2  s-1. 

Dh2  0,m 

diffusivity  of  water  in  the  membrane,  m2 

ex,  ey,  ez 

coordinate  vectors 

Estack 

stack  voltage,  V 

F 

Faraday’s  constant,  A  s  mol-1 

i,  i 

current  density,  A  m  2 

j 

building  block 

k 

thermal  conductivity,  W  m-1  K-1 

Mm 

equivalent  weight  of  the  dry  membrane, 

n 

the  number  of  cells 

no 

electroosmotic  drag  coefficient 

N> 

molar  flux  of  species  1,  mol  m-2  s-1 

P 

pressure,  Pa 

S 

source  term 

T 

temperature,  K 

u,v,v,U 

velocities,  m  s-1 

x,y,z 

coordinate,  m 

Greek 

K 

permeability,  m2 

A 

water  content 

P 

dynamic  viscosity,  kg  m-1  s-1 

P 

density,  kg  m  3 

<7m 

protonic  conductivity,  S  m-1 

0's 

electric  conductivity,  S  m  1 

<P 

potential,  V 

stream  function 

Superscripts 
cool  cooling 

eff  effective 

(1)  liquid 

Subscripts 
H2  hydrogen 

H20  water 

i  species  i 

m  membrane 


mom  momentum 

N2  nitrogen 

02  oxygen 

pot  potential 

s  solid  phase 

temp  temperature 
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